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Abstract. The solar neutrino experiment Borexino, which is located in the Gran Sasso 
underground laboratories, is in a unique position to study muon-induced backgrounds in 
an organic liquid scintillator. In this study, a large sample of cosmic muons is identified 
and tracked by a muon veto detector external to the liquid scintillator, and by the specific 
light patterns observed when muons cross the scintillator volume. The yield of muon-induced 
neutrons is found to be Y n = (3.10±0.11) • 10~ 4 n/(fi- (g/cm 2 )). The distance profile between 
the parent muon track and the neutron capture point has the average value A = (81.5 ± 
2.7) cm. Additionally the yields of a number of cosmogenic radioisotopes are measured for 
12 N, 12 B, 8 He, 9 C, 9 Li, 8 B, 6 He, 8 Li, n Be, 10 C and n C. All results are compared with Monte 
Carlo simulation predictions using the Fluka and Geant4 packages. General agreement 
between data and simulation is observed for the cosmogenic production yields with a few 
exceptions, the most prominent case being 11 C yield for which both codes return about 50% 
lower values. The predicted [i-n distance profile and the neutron multiplicity distribution are 
found to be overall consistent with data. 

Keywords: Borexino, Muon, Cosmic, Cosmogenic, Neutron 
1 Introduction 

The Borexino experiment is a 278 1 liquid-scintillator detector designed for real-time mea- 
surements of low energy (<20MeV) neutrinos. The primary goal of the experiment is the 
spectroscopy of solar neutrinos. In this respect Borexino has performed a precision measure- 
ment of the Be neutrino line [1, 2], has lowered the threshold for the real-time detection of 
the 8 B neutrino spectrum to 3MeV [3], and has directly observed the neutrinos of the pep 
line, at the same time placing the most stringent limit on the CNO neutrino flux [4] . Borex- 
ino is also very competitive in the detection of anti- neutrinos having reported a first 
observation of geo- neutrinos in 2010 [5], followed by a recent new measurement [6]. Finally, 
the experiment is sensitive to neutrino signals from a galactic core-collapse supernova [7, 8]. 

The extremely low background in the scintillator target is essential to the success of 
Borexino. While careful pre-selection of detector materials and extensive purification of the 
organic scintillator is necessary, shielding from external, and especially cosmic, radiation is 
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of comparable importance. The detector is located deep underground (3800 meters of water- 
equivalent, m w.e.) in the Hall C of the Laboratori Nazionali del Gran Sasso (LNGS, Italy), 
where the cosmic muon flux is suppressed by about six orders of magnitude. In spite of this 
large attenuation factor, residual muons constitute an important source of background for 
neutrino detection. For example, they produce neutrons or radioactive isotopes by spallation 
reactions in target materials, e.g. 12 C, and these produce signals which mimic the observation 
of a reaction of interest. 

An understanding and mitigation of muon-induced backgrounds are of great relevance 
to all investigations of rare processes. In the majority of underground rare-event experiments, 
muons and their spallation products constitute a severe source of background. For example, 
in direct dark matter searches, neutron interactions feature a signature very similar to those 
induced by WIMPs, and careful shielding must be employed. In 0vf3f3 experiments, the /3- 
decays of long-lived radioisotopes produced in-situ can be significant background components. 
It is expected that cosmogenic backgrounds will be even more important in the next gener- 
ation of low-background experiments, as the detector sizes and sensitivities are increasing. 
Thus more sophisticated muon vetoes as well as extensive shielding will be necessary. 

This paper presents the results of a detailed investigation of muons, and especially their 
spallation products, which were detected by Borexino. Due to its simple geometry, large 
mass, and excellent event reconstruction capability, Borexino offers the unique possibility for a 
precise study of cosmogenic production inside a large, uniform volume of low-Z material. The 
paper is structured as follows. After introducing the detector layout and the general cosmic 
background conditions found at the LNGS facility (section 2), we review the results of the 
cosmic muon flux and present the reconstructed angular distribution of the muons (section 3). 
We continue with a detailed study of the neutron production rate and multiplicity in liquid 
scintillator (section 4). We also present the lateral distance profile of neutron captures with 
respect to the parent muon track. This is of special interest for dark matter experiments 
relying on low-Z materials for neutron shielding. Moreover, we investigate the production 
of a selection of cosmogenic radioisotopes in the scintillator volume (section 5). Finally, 
we perform a detailed comparison of our experimental results with the rates and profiles 
predicted by the commonly used Fluka and Geant4 simulation codes (section 6). This 
validity check is of considerable importance as Monte Carlo simulations represent virtually the 
only means to transfer our findings to the various detector geometries realized in other low- 
background experiments. We also show the production yields determined in the KamLAND 
experiment [9] which features a setup largely comparable to Borexino. Section 7 summarizes 
our main results. 

2 The Borexino Detector at the LNGS 

The Borexino detector [10] consists of a spherical inner detector (ID) containing a liquid 
scintillator target and a surrounding outer detector (OD) consisting of a large water tank. 
This tank acts both as passive shield and as an active muon veto. The general layout is 
presented in figure 1. The central, active scintillator consists of pseudocumene (PC, 1,2,4 
trimethylbenzene) , doped with 1.5g/liter of PPO (2,5-diphenyloxazole, a fluorescent dye). 
The nominal target mass is 278 1. The scintillator is contained in a thin (125 pm) nylon 
vessel of 4.25 m radius and is shielded by two concentric inactive PC buffers (323 1 and 567 1) 
doped with few g/1 of a scintillation light quencher (dimethylphthalate). The two PC buffers 
are separated by a second thin nylon membrane to prevent diffusion of radon towards the 
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Figure 1. Sketch of the Borexino detector. 



scintillator. The scintillator and buffers are contained in a Stainless Steel Sphere (SSS) with 
a diameter of 13.7 m. The SSS is enclosed in a 18.0 m diameter, 16.9 m high domed Water 
Tank (WT), containing 2100 1 of ultra-pure water. The scintillation light is detected via 
2212 8"-photomultiplier tubes (PMTs) uniformly distributed on the inner surface of the SSS. 
Additional 208 8" PMTs instrument the WT and detect the Cherenkov light radiated by 
muons in the water shield [11]. 

In Borexino, low energy neutrinos (u) of all flavors are detected by means of their elastic 
scattering off electrons or, in case of electron anti-neutrinos (P e ), via the inverse f3 decay on 
free protons. The electron (positron) recoil energy is converted into scintillation light which 
is then collected by the ID PMTs. Borexino is sensitive to neutrinos of at least 100 keV in 
energy, while the inverse f3 decay induced by anti-neutrinos requires a minimum neutrino 
energy of 1.8 MeV. While cosmic muons crossing the ID deposit much greater energies and 
create substantially more light, cosmogenic neutrons and radioisotopes induce scintillation 
signals on a scale similar to neutrino interactions. 

3 Cosmic Muons 

The primary cosmic muon flux arriving at Earth's surface (6.5-10 5 /j,/(m 2 -h)) is strongly 
attenuated when penetrating the mountain above the detector by about a factor 10 6 . The 
rock shielding is equivalent to some 3800 m of water. Thus the mean energy of the muons at 
LNGS site is about 280 GeV, compared to about 1 GeV at the surface, since the lower energy 
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Figure 2. Angular distribution of the muons crossing the Borexino IV. The side plots show the 
projections to astronomic azimuth and zenith angles. 



muons incident at the surface are absorbed, and the spectrum steeply falls as a function of 
energy [12]. 

We consider in this article only muons generating signals in both the ID and the OD for 
which track reconstruction is also performed. Muons identification occurs with three different 
methods [11]. The first two are based on the detection of the Cherenkov light produced in 
the water. In the first method, the light triggers the OD sub-system (muon trigger flag, 
MTF). In the second method, a cluster of OD PMT pulses is identified, correlated in space 
and time (muon cluster flag, MCF). The third method relies on pulse shape identification 
of muon tracks among the point-like scintillation events detected by the ID (inner detector 
flag, IDF). The detection efficiencies are 0.9925(2), 0.9928(2) and 0.9890(1) respectively. The 
cosmic muon interaction rate in Borexino was found in [13] to be (4310 ± 2 sta t ± 10 sys t) d _1 
and corresponds to a flux of (3.41 ± 0.01) • 10~ 4 m~ 2 s _1 as measured in Hall C of the LNGS 
laboratory. 

The Borexino ID features a uniform acceptance for incident cosmic muons which is 
independent of the arrival directions, thanks to its spherical symmetry. The observed angular 
distribution is shown in figure 2 as a function of the astronomic azimuth (</>) and zenith (8) 
angles in a two-dimensional contour plot 1 . In addition, the one-dimensional projections on 
the 9 and 0-planes are also presented. All three distributions reflect the influence of the local 
mountain topology: The differences in the thickness of the overlaying rock are imprinted as 
angle-dependent variations in the residual muon flux. Note that due to its uniform detection 

A table listing the numerical values of the distribution has been included in the supplementary materials. 
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efficiency, Borexino is in a unique position to map the muon distribution at depth close to 
the horizon without angular distortions. 

To obtain these plots, the ID tracking algorithm described in [11] was employed. These 
plots consider only muons crossing the IV because their tracks offer the best angular resolu- 
tion. For this selection, the (redundant) requirements of a reconstructed energy deposition 
of at least 300 MeV in the IV and an impact parameter of less than 4.25 m were applied. The 
remaining sample consists of 1 221470 individual muon tracks. 

4 Cosmogenic Neutrons 

Cosmic muons crossing the detector produce fast neutrons through different spallation pro- 
cesses on carbon nuclei. The neutron velocities are slowed in the scintillator by collisions 
with hydrogen or carbon nuclei to thermal sub-eV energies in a few scatterings. This process 
occurs within a few tens of ns. Consequently, signals from ionizations due to the recoiling 
nuclei cannot be disentangled from the much higher light emission of an incident muon. For 
the same reason, fast neutron captures are also not visible. About 1% of the neutrons are 
captured during a fast capture process, as determined by our Monte Carlo simulations. We 
do not correct our results to account for this effect in this paper, and hereafter we refer to 
all neutron captures as due to slower thermal capture. The mean capture time of a ther- 
mal neutron in the liquid scintillator is ~250us, and the subsequent gamma-ray emission is 
distinctly visible. The energy emitted in gamma rays is 2.2 MeV if the neutron is captured 
on hydrogen and 4.9 MeV if captured on carbon. Based on the elemental composition of the 
scintillator and the relative capture cross-sections, about 99% of all thermal neutron captures 
are expected on hydrogen [14]. The ability of Borexino to detect cosmogenic neutrons was 
described in detail in [11]. 

4.1 Neutron detection with the main electronics 

The data acquisition system issues a dedicated neutron acquisition gate of 1.6 ms length after 
a muon crosses the SSS. The captured PMT pulses are due to neutron capture-gammas and 
accidental 14 C decays. The latter is an intrinsic contaminant of the scintillator. Particular 
care must be taken when the muon crosses the scintillator because of the large amount of 
visible light which is created. For these events, the baseline of the front-end electronics takes 
up to 30 us to stabilize. Furthermore as a result of the intense PMT illumination, the front 
of the acquisition gate is highly populated by noise pulses creating a variable baseline on 
which the other pulses are superimposed. Thus the active time window for the analysis 
is set after the baseline stabilizes (30 us). An ad- hoc algorithm detects physical events by 
identifying clusters of time-correlated PMT pulses on top of the baseline with high efficiency. 
An energy threshold of 1.3 MeV is imposed via a variable selection cut based on detector 
saturation. The overall detection efficiency above threshold was determined using neutron 
samples which were selected as described below. The result is £det = (91.7±1.7 sta t±0.9 sys t) % 
for cosmogenic neutrons captured later than 30 us after their parent muon is observed. The 
error incorporates the uncertainty in the threshold definition of the energy. 

The data sample used for the analysis contains 559 live days, covering the time period 
from January 6, 2008, until February 2, 2010. The respective average scintillator volume 
contained in the nylon vessel was measured to be (306.9 ± 2.9) m 3 , or (22.8 ± 0.2)% of 
the SSS volume. This translates into an active mass of (270.2 ± 2.6) t. The volume in 
which neutron captures are detected is not identical with that defined by the physical vessel 



- 5 - 



boundaries. For example, gamma-rays generated in the buffer close to the nylon vessel may 
cross into the scintillator volume and create a sizeable light output, while gamma-rays from 
neutron captures inside the scintillator may escape into the buffer undetected. This effect 
was quantified using a Monte Carlo simulation of neutron capture gammas in both buffer 
and scintillator. The simulation included the transport of scintillation light and electronic 
effects. Of all 7's generated inside the SSS, the fraction of events which deposit more than 
1.3 MeV energy in the IV, £ vo i, is (23.0 ± 0.3) % for captures on hydrogen, and (24.2 ± 0.2) % 
for captures on carbon. Both the physical volume evaluation and the Monte Carlo simulation 
take into account variations in actual vessel shape and size over the data collection period. 

The neutron capture time and the purity of the sample were studied based on the time 
difference (At) between the occurrence of candidate clusters and their parent muons. The At 
distribution is fit by the sum of an exponential decay and a flat component for uncorrelated 
events. The resulting neutron capture time is r n = (259.7 ± 1.3 s tat ± 2.0 sys t)us, and is in 
agreement with our previous measurement [11]. Based on this value, the fraction of neutrons 
captured later than 30 us after the parent muon is £t = (89.1 ± 0.8) %. The contamination 
by uncorrelated events is found to be (0.5 ± 0.2) %. 

4.2 Neutron detection with waveform digitizers 

The Borexino detector is equipped with auxiliary DAQ systems based on fast waveform digi- 
tizers. Two of these systems (hereafter SYS1 and SYS2) were used in this analysis to evaluate 
the neutron detection efficiency £det of the main data acquisition system. Furthermore, these 
two systems provide a cross check on the neutron yield measurement (section 4.3). 

SYS1 is a single channel 500 MHz, 8 bit digitizer (Acqiris DP235) which records the 
cumulative analog output of all ID PMTs. It is triggered by the MTF condition of the outer 
detector, and collects data for about 1.6 ms. A cluster-finding algorithm identifies gamma-ray 
capture signals between 30 and 1590 us after muon detection. An energy threshold of 1.3 MeV 
is applied to reject noise pulses. SYS1 collected data between April 2008 and November 2009. 
The fit to the At distribution returns a capture time of (261 ± l s t a t ± Isyst) us. An additional 
flat component which would be allowed by the fit is consistent with zero at the lcr confidence 
level. 

SYS2 is an auxiliary hardware architecture, based upon 96 waveform digitizers (CAEN 
v896: 400 MHz, 8 bit). Each channel receives the analog sum of 16 or 24 PMTs. The system 
operates independently from the main DAQ. A separate trigger is implemented by an FPGA 
unit (CAEN vl495). Data used for this analysis were collected between December 2009 and 
March 2012. A neutron detection efficiency of essentially 100% is reached for an energy 
threshold above 1 MeV. For SYS2, the fit to the At distribution returns a capture time of 
(258.7±0.8 s tat±2.0 sys t) ps. A residual of uncorrelated background on the level of (0.5±0.1) % 
was determined by a fit and when using a delayed time window. 

The values of the capture time obtained by SYS1 and SYS2 are in good agreement with 
the values using the main electronics (section 4.1). 

4.3 Neutron production rate and multiplicity 

The data set acquired by the main DAQ contains a sample of = 2 384 738 muons and N n = 
111 145 neutrons. The neutron capture rate for the efficiencies described above, is determined 
to Rn = (90.2 ± 2.0 s tat ± 2.4 sys t) (d 100t) -1 after scaling. The systematic uncertainty of this 
measurement is dominated by the neutron capture time and the average scintillator volume 
contained inside the IV. 
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Figure 3. The muon-neutron distance distribution observed in Borexino: black crosses represent 
the data points for the standard neutron hit multiplicity cut. The shaded-grey area indicates the 
systematic uncertainty. The fit of the toy Monte Carlo is indicated by the solid red line. The dashed 
lines correspond to two exponential components, each featuring a decay length A and a relative fraction 
/. The muon resolution parameters \i and a are left free in the fit procedure (see sect. 4.4 and [11] 
for details). The table lists the best-fit results with statistical and systematic uncertainties. The fit 
returns x 2 / n df = 57/54. 



The rate of muons which produce neutrons that eventually are captured inside the 
IV, is Rn = (67.5 ± 0.4 sta t ± 0.2 syst ) d^ 1 (~1.5% of muons crossing the ID). The detector- 
specific ratio R n / Rn corresponds to an average neutron multiplicity within the IV volume of 
M = (3.61 ±0.08 stat ±0.07 syst ) n/ fi. The distribution of the multiplicities of detected neutron 
captures is shown in figure 13 where it is compared to Monte Carlo predictions. The neutron 
yield per unit length of muon track in the target medium is 

v J_ _i 1 

n AT ' /? av S 

^ V C/u Pscint £det " e t ' £yol 

= (3.10 ± 0.07 sta t ± 0.08 syst ) • 10~ 4 n/Ou • (g/cm 2 )) (4.1) 

where £^ g = 4/3i?sss is the average muon path through the SSS with i?sss = (6.821 ± 
0.005) m and 

Pscint — 0.88 g/cm^ is the scintillator density. We chose to consider the muon 
path through the SSS and not just the IV and to include the ratio of the two volumes in e vo i 
in order to correctly account for the effective neutron detection volume. This is discussed 
in section 4.1. The statistical uncertainty associated with this result has been assessed by 
a toy Monte Carlo code which simulates neutron production by muons in order to observe 
the size of the fluctuations. For varying values of the muon rate, the statistical uncertainty 
is approximately seven times the square root of the number of neutron captures divided by 
the live time. This results in a statistical uncertainty of 3%. 

As a consistency check, the neutron yield has been also determined based on the wave- 
form digitizers: V n SYS1 = 3.19±0.08 stat i^ syst for SYSl and V n SYS2 = 2.87±0.07 stat ±0.15 syst 
for SYS2 in units of 10 _4 n/(^ • (g/cm 2 )). This is in reasonable agreement with the value 
given in equation (4.1). 

4.4 Neutron lateral distance 

Borexino can spatially reconstruct the emission point of the neutron capture gamma-ray as 
well as the track of its parent muon. This information may be used to compute the shortest 
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(perpendicular) distance between the neutron vertex and the parent muon track (hereafter 
denoted as the neutron lateral distance). A toy Monte Carlo was fit to the data in order to 
separate the distance travelled by the neutrons from resolution effects of the detector (see 
below) . 

Muon track and neutron verticies were reconstructed based on the Borexino main DAQ 
in order to obtain the lateral distance distribution from the data. Only muons and neutrons 
having a radial distance less than 4 m from the detector center were selected in order to en- 
sure a well-defined geometry for comparison to the Monte Carlo. Moreover, the samples were 
cleaned to remove tracks and vertices of inferior reconstruction quality. In case of the muons, 
only events which feature spatially compatible tracks in both sub-detectors are used. The 
neutron selection cuts are much more restrictive. Due to electronic effects and PMT after- 
pulses which were present after very luminous muons, the spatial reconstruction of subsequent 
neutron events can be severely compromised. This results in systematic shifts increasing or 
decreasing the distance between the neutron capture point and the parent muon track [11]. 
Systematic studies demonstrated that the majority of these effects (e.g. afterpulses) subside 
in the first 200 us after a muon event. At later times, a fraction of the electronic channels 
might be affected by buffer overflow, which leads to asymmetric PMT hit patterns. Therefore, 
only neutron captures with a minimum time delay of 200 us are compared to a parent muon 
with the further restriction that the neutron event was composed of at least 100 individual 
PMT signals in order to ensure an unbiased vertex reconstruction (iVhits > 100). The latter 
condition is easily met by neutrons produced by minimum-ionizing muons, but depletes the 
sample of useful neutron captures in case of very luminous muon events. Such muons are 
expected to create extensive hadronic showers, and there is a risk that the hit multiplicity 
cut for the neutrons introduces a bias to the selected sample which prefferentially suppresses 
neutrons at large distances from their parent tracks. Finally, we limit the visible energy 
window to .E v i s 6 [0.9; 4.8] MeV in order to select only neutron captures on hydrogen and 
carbon, while removing a minor contamination from short-lived cosmogenic isotopes. The 
combination of cuts reduces the remaining sample to ~20 % of the original neutrons. 

The resulting lateral distance distribution is shown in figure 3. The grey shaded area 
corresponds to the systematic uncertainty introduced by the cut N^ its > 100, and was ob- 
tained by varying the minimum iVhits condition for neutron selection from to 200. Due to 
the broad energy spectrum of the spallation neutrons and the corresponding distribution of 
the neutron mean free paths, it is found that two exponential components (A s hort, and Ai ong ) 
are satisfactory for an empirical representation of the distribution. The fit function shown 
in figure 3 was obtained by a toy Monte Carlo simulation. Apart from the exponential com- 
ponents, the fit takes into account the muon and neutron spatial resolutions, which includes 
the average displacement of the neutrons during thermalization and the finite propagation 
distance of the capture gamma in scintillator (~20cm). The geometric impact of the ap- 
plied radial cuts described above are included. The muon lateral resolution is described by 
a Gaussian smearing a with a constant radial offset \x. These are free parameters in the fit. 
Conversely, the neutron vertex resolution is set to a fixed value of 23cm (see [11] for details). 

The fit returns a short component A s hort = (61.2 ± 0.6 sta t ± 2.6 syst ) cm which is in 
agreement with earlier LVD results [15]. The long component is found to be Ai ong = (147 ± 
3 sta t i 12 sys t) cm. Systematic uncertainties for the parameters were determined by multiple 
repetitions of the fit while varying the minimum A^its condition for the neutrons. Based 
on the relative weights of the two effective components, an average lateral distance of A = 
(81.5 ± 2.7) cm was determined. 
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Table 1. List of cosmogenic isotopes expected to be produced by muons in organic scintillators in 
measurable rates. 

5 Cosmogenic Radioisotopes 

In addition to neutrons, radioactive isotopes are produced in muon- induced spallation pro- 
cesses on the target nuclei. A list of the relevant cosmogenic isotopes with their properties 
and sorted by increasing lifetime can be found in table 1. 

Candidate events are selected via two observables: the visible energy, E, and the time 
difference, At, with respect to a parent muon. The distributions of the two observables 
are fit simultaneously in an unbinned likelihood fit. However, matching to a parent muon 
is not unique in general, as many muons can be present within the selected analysis time 
gate t g . This results in multiple values of At for a given candidate, with only one muon 
physically correlated with the neutron decay. This effect is most prominent for the analyses 
on cosmogenic isotopes with lifetimes on the order of seconds or longer since in average cosmic 
muons cross the ID every 20 s. The distribution in At is fit with the function: 

AT* At jv* N* 
F(At) = Y,—e~^ + — + — (5-1) 

i T i tg tg 

The number of decays in each isotope profile is Nj with r% its lifetime. Flat contributions, 
Nt and iV* m , account for uncorrelated background events and for physically uncorrelated 
matches, respectively. The latter is a property of the selected data set and calculated in- 
dependently. The fit function is valid for time scales much shorter than the average run 
duration of a data set (~6h). This is valid for all cosmogenic isotopes with the exception 
of n C (table 1). As will be presented in section 5.7, a distortion of the time profile due to 
run-boundary effects can be avoided by a time cut of candidate events which occur close to 
run-start. The spectral shapes of the respective isotopes are generated with the Geant4 
based Borexino Monte Carlo code. The simulation reproduces the full detector response, 
yielding about 500photoelectrons/MeV of deposited energy in j3 and 7 decays, and spectral 
fits can be performed directly on the number distribution of photoelectrons from candidate 
events. For easier reference, this conversion factor will be used to refer to the energy selection 
cuts in the [MeV] units in the following analyses. Based on the visible energy E of an event, 
the spectral fit function G(E) is given by : 

G(E) = £ Nf m (E) + N b E g b (E) (5.2) 

i 

The spectral shapes of the analyzed isotopes are denoted by g%(E) with the respective number 
of decays Nf. The uncorrelated background is addressed with the spectral shape gb(E) and 
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Figure 4. Simultaneous fit of the cosmogenic isotopes 12 N and 12 B in visible energy deposition 
(top panel) and decay time relative to preceding muons (bottom panel), including the isotopes 8 He 
(blue), 9 C, 9 Li, 8 B and 8 Li as contaminants. The fit returns only upper limits for the isotopes 12 N, 
8 He, 9 C and 8 Li, and they cannot be seen in the graph. The goodness of the simultaneous fit is 
X 2 /ndf= 348/236. 



is the number of entries. The energy distribution is generated from events which occur 
within a time interval, ig, relative to preceding muons. To enhance the signal to noise ratio, 
tE is chosen to be in the order of the lifetimes of the respective cosmogenic radionuclides. 
The number of isotope decays in the time and energy fit (Nf and Nf), and the background 
events (N£ and iVf), are related to the known selection cut efficiencies in time and energy. 

Most cosmogenic isotopes are expected to be produced at a very low rate. To reduce 
accidental coincidences, muons are removed from the sample of candidate isotope events by 
the application of the MTF and IDF muon identification methods (section 3). The efficiency 
and small distortions due to the application of IDF are included in the Monte Carlo generated 
spectral shapes. MTF efficiency is accounted for a posteriori. Unless stated otherwise, the 
positions of candidate events are also required to lie within a fiducial volume (FV). This is 
defined by a sphere of a 3m radius which corresponds to a 99.6 1 mass of liquid scintillator. 
The systematic uncertainty in the reconstructed volume for the decay energies of interest is 
estimated to ± 3.8 % and contributes to the uncertainty of all measured yields. Except where 
differently noted, the analyses on cosmogenic radioisotopes are based on the same data set 
used for the neutron yield analysis (section 4.1). 

5.1 12 N and 12 B 

Candidate events for the decays of 12 N (/3 + -emitter, r = 15.9 ms, Q = 17.3 MeV) and 
12 B (/^"-emitter, r = 29.1ms, Q = 13.4 MeV) are selected within an energy range E 6 
[3.6, 18] MeV and a time gate At £ t g = [2 ms, 10 s] with respect to a preceding muon event. 
The energy distribution is constructed from events with At G tE = [2,60]ms. This in- 
creases the signal-to-background ratio. After each muon, a 2 ms veto is applied to avoid 
muon-induced secondaries (mainly neutrons) which leads to a negligible dead time. Fur- 
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Lateral Distance to Muon Track (cm) 

Figure 5. Comparison of measured and simulated lateral production profiles for cosmogenic 12 B 
candidates inside the FV with respect to the parent muon track. For improved spatial resolution we 
select only tracks of muons which cross the IV at an impact parameter of less than 4 m. 



thermore, decays of the cosmogenic isotopes 8 He, 9 C, 9 Li, 8 B and 8 Li are considered as 
contaminations and are fit alongside 12 N and 12 B. The upper limit of the time gate t g (i.e. 
10s) is driven by the lifetimes of the isotopes 8 B (r = 1.11s) and 8 Li (r = 1.21s). In 
addition, a fraction of (86.2 ± 0.2) % of all n Be decays is expected within the selected 
energy range. However, due to its low production rate and long lifetime (r = 19.9 s), 
this contribution is estimated to be less than 1% in the time profile, and negligible in 
the energy distribution (section 5.5). The background spectral shape is built from events 
with At > 10 s to avoid accidental coincidences of short-lived cosmogenic isotopes. Fig- 
ure 4 shows the simultaneous fit in energy and time. The efficiency of the energy cut 
is evaluated via the Borexino Monte Carlo simulation to be e( 12 N) = (79.3 ± 0.4) % and 
e( 12 B) = (84.0 ± 0.3) %. The uncertainties are due to the detector energy resolution. The 
simultaneous fit yields a rate of R( 12 N) < 0.03 (d 100t) _1 at a 3a confidence level, and 
i2( 12 B) = (1.62 ± 0.07 stat ± 0.06 syst ) (d lOOt)" 1 . These translate to production yields of 
Y( 12 N) < l.l-10- 7 /(Mg/cm 2 )) andV( 12 B) = (55.6 ± 2.5 stat ± 2.1 syst ) • 10" 7 /(p ■ (g/cm 2 )). 

12 B at a level of (94.6 ±0.3) % of all events is clearly the dominating cosmogenic isotope 
within the selected time window (At E [2,60] ms). Figure 5 shows the lateral production 
profile of this isotope with respect to the reconstructed track of the parent muon. 

5.2 8 He and 9 Li 

Both /^--emitters 8 He (r = 171.7 ms , Q = 10.7 MeV) and 9 Li (r = 257.2 ms, Q = 13.6 MeV) 
exhibit daughter nuclei with neutron-unstable excited states. With a 16% branching ratio, 
the /?-decay of 8 He populates such a state in 8 Li. For 9 Li, the branching ratio to a neutron- 
unstable state in 9 Be is 51%. The subsequently emitted neutron is captured mainly on 
hydrogen with a mean capture time of (259.7±3.3) us (section 4.1), emitting a characteristic 
2.2 MeV gamma-ray. The triple-coincidence of a muon, a /3-emission, and a delayed neutron 
capture provides a very clean signature, which allows analysis of events within the entire 
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Figure 6. Simultaneous fit of the cosmogenic isotopes 8 He and 9 Li in visible energy deposition (top 
panel) and decay time relative to preceding muons (bottom panel). The fit returns only an upper 
limit for the isotope 8 He and its line cannot be seen in the graph. The goodness of the simultaneous 
fit is x 2 / ndf = 71/98. 



mass of the inner vessel. Thus, an enlarged data set of 1366 live days taken between January 
6, 2008, and August 31, 2012, is used for the analysis with a mean inner vessel mass of 
(268.2 ± 2.8) t. The requirements for candidate events for the /3-emissions are At G [2 ms, 2 s] 
and E G [0.8, 14] MeV. The energy distribution is taken from events with At G [2ms,ls]. 
Subsequent neutron capture candidates are selected by E G [1.7, 2.6] MeV. These are required 
to occur within lm distance and a maximum time delay of 1.3 ms to a /3-like event. The 
uncorrelated background spectrum is derived from events with At > Is. Figure 6 shows the 
simultaneous fit in energy and time of 8 He and 9 Li. The j3n selection cut efficiency has been 
evaluated to be e{j3n) = (79.3 ±0.4) % via the Borexino Monte Carlo simulation. The energy 
cut efficiencies are estimated to e( 8 He) = (99.49 ±0.05) % and e( 9 Li) = (96.99 ±0.11) %. The 
fit returns only an upper limit for 8 He. Isotope production rates of i?( 8 He) < 0.042 (d lOOt)™ 1 
at a 3a confidence level, and R( 9 U) = (0.083 ± 0.009 sta t ± 0.001 syst ) (d lOOt)" 1 are observed. 
The corresponding yields are y( 8 He) < 1.5 • 10~ 7 /(/i • (g/cm 2 )) and y( 9 Li) = (2.9 ± 0.3) • 
10- 7 /(/"-(g/cm 2 )). 

5.3 8 B, 6 He and 8 Li 

The cosmogenic isotopes 8 B (/3 + -emitter, r = 1.11s, Q = 18.0 MeV), 6 He (/3 _ -emitter, 
r = 1.16 s, Q = 3.51 MeV) and 8 Li (/^--emitter, r = 1.21s, Q = 16.0 MeV) feature similar 
lifetimes. However, the significantly lower Q-value of 6 He enables a partial disentanglement 
of these radionuclides via cuts in visible energy and time. We separate the energy range in 
two regimes, denoted as ER1 (E £ [2, 3.2] MeV) and ER2 (E G [5,16] MeV), respectively. 
Regime ER1 comprises decays of all three isotopes, whereas ER2 includes only 8 B and 8 Li. 
The two energy intervals are fit simultaneously with their respective time profiles in a single, 
un-binned maximum likelihood fit. The spectral shape of uncorrelated background is derived 
from events with At > 140 s. Table 2 summarizes the energy selection efficiencies and the 
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Energy regime 


j. ime Lxate tg 


lime Lraie le 


Cosmogenic 


Lifetime 


Energy Cut 




M 


[s] 


Isotope 


M 


Efficiency [%] 








8 B 


1.11 


4.0 ±0.3 


ERl 


[1, 140] 


[1,2] 


6 He 


1.16 


16.8 ±0.3 


E e [2,3.2]MeV 






8 Li 


1.21 


8.1 ±0.2 








io c 


27.8 


77.1 ± 0.2 


ER2 


[1,10] 


[1,3] 


8 B 


1.11 


81.6 ±0.4 


E e [5, 16] MeV 






8 Li 


1.21 


67.5 ±0.4 



Table 2. Selection cuts of the candidate events in energy (E) and time (t g , £e) for the two regimes ERl 
and ER2 for isotopes 8 B, 6 He and 8 Li. In addition, the expected cosmogenic isotopes are given with 
their lifetimes and energy cut efficiencies. For each energy regime, the time profile is constructed from 
events within the time interval t g relative to preceding muons. To enhance the signal-to-background 
ratio, the energy distribution is based on events in the time interval tE- The distributions and the 
result of a simultaneous fit in time and energy to both regimes is shown in figure 7. 

chosen time gates for the time and energy distributions in these two energy regimes. Due 
to the lower energy threshold of ERl, an additional contribution of the cosmogenic isotope 
10 C (/3 + -emitter, r = 27.8 s, Q = 3.65 MeV) is included as a free parameter in the fit. 
To avoid contaminations of short-living cosmogenic isotopes, a 1 s veto after each muon is 
applied for both regimes, inducing a dead time of 3.6 %. The result of the simultaneous 
fit is shown in figure 7 for the energy regimes ERl and ER2, respectively. The isotope 
production rates are found to be R( S B) = (0.41 ± 0.16 sta t ± 0.03 syst ) (dlOOt)" 1 , i2( 6 He) = 
(1.11 ± 0.45 stat ± 0.04 syst ) (d lOOt)" 1 and R{ s Li) = (0.21 ± 0.19 stat ± 0.02 syst ) (d lOOt)" 1 . The 
corresponding yields are Y( 8 B) = (1.4±0.6 sta t±0.1 syst )-10- 6 /V (g/cm 2 )), Y( 6 He) = (3.80± 
1.53 stat ±0.14 syst )-10- 6 /(^-(g/cm 2 )) and Y( 8 Li) = (7.1±6.6 stat ±0.7 syst )-10" 7 /(/i- (g/cm 2 )). 

5.4 9 C 

Candidate events for the decay of 9 C (/3 + -emitter, r = 182.5 ms, Q = 16.5 MeV) are selected 
by E £ [5, 18] MeV and At £ [250 ms, 10.25 s]. The lower energy threshold avoids decays 
of 6 He in the data set. Events occurring for At S [250,600] ms are employed to obtain the 
energy distribution. After each muon, a 250 ms veto is applied, rejecting contributions from 
shorter-lived cosmogenic radionuclides and reducing the live time of the data set by 0.9%. 
The isotopes 8 He, 9 Li, 8 B and 8 Li are treated as contaminants. Events for At > 10.25 s are 
used to build the spectral shape of the uncorrelated background. The energy distribution 
for the 8 B, 6 He and 8 Li analysis (ER2 with E £ [5, 16] MeV, t E € [1, 3] s) is used to confine 
the rates of 8 B and 8 Li in the simultaneous fit as additional complementary information 
(section 5.3). 

The result of the best fit for the time profile and both energy distributions is shown in 
figure 8. A fraction of e( 9 C) = (73.4 ± 0.4) % of all 9 C candidates is expected within the 
energy range of the 9 C candidates. The rate and yield of 9 C are determined to an upper 
limit of R( 9 C) < 0.47 (d lOOt)" 1 , and Y( 9 C) < 1.6 • 10" 6 /(// • (g/cm 2 )) at 3a confidence level 
after correcting for efficiencies. 

5.5 n Be 

Events with E £ [5, 12] MeV are chosen in order to determine production rate for n Be which 
is a f3~ -emitter , with r = 19.9 s, Q = 11.5 MeV. The lower energy boundary is needed in order 
to reject 10 C and 11 C decays as well as the external (uncorrelated) 7-background from 208 T1. 
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Figure 7. Simultaneous fit of the cosmogenic isotopes 8 B, 6 He and 8 Li in visible energy deposition 
(first and third panels) and decay time relative to preceding muons (second and fourth panels) for 
the energy regimes ER1 (E G [2, 3.2] MeV) and ER2 (E G [5, 16] MeV). The isotope 10 C is included 
as contaminant. Some of the isotopes cannot be seen in the graph because the fit returns a very low 
value for their rates. The goodness of the simultaneous fit is % 2 /ndf = 457/499. 



The time profile is constructed for At G [10, 210] s with respect to a preceding muon whose 
track is reconstructed within 1.5 m from the n Be candidate decay. The energy distribution 
is composed of events with At £ [10, 40] s in order to increase the signal-to-background ratio 
after applying the same muon track cut. A 10 s veto after each muon rejects the shorter-lived 
cosmogenic radionuclides, which decreases the live time of the data set by 28.4 % and leaves 
11 Be as the only cosmogenic isotope. The background spectral shape is derived from events 
which satisfy the muon track cut and occur later than 210 s after the muon. The best fit 
results are shown in figure 9. The selection efficiency of the muon track cut is estimated to 
(63.3 ± 2.5) % using the time and lateral distribution of cosmogenic 12 B to preceding muon 
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Figure 8. Simultaneous fit of the cosmogenic isotope 9 C in visible energy deposition (top panel) and 
decay time relative to preceding muons (bottom panel). The fit returns only an upper limit for the 
isotope 9 C and some isotopes which cannot be seen in the graph. The goodness of the simultaneous 
fit is x 2 / ndf = 218/268. 



tracks. The fraction of decays in the selected visible energy range is calculated to to be 
e( 11 Be) = (69.3 ± 0.6) %. The isotopic production rate and yield of n Be are measured to be 
R( n Be) < 0.20 (d lOOt)" 1 and Y( n Be) < 7.0 • 10" 7 /(p • (g/cm 2 )) at the 3ct confidence level. 

5.6 10 C 

The production of 10 C (/3+-emitter, r = 27.8 s, Q = 3.65 MeV) in muon-induced spalla- 
tion processes on 12 C is usually accompanied by the emission of at least one free neutron. 
These neutrons are eventually captured on hydrogen or carbon (section 5.2). The acciden- 
tal background is significantly reduced after requiring a three-fold coincidence between a 
muon, at least one subsequent neutron capture, and a 10 C decay candidate. Neutron cap- 
tures with a minimum energy of 1.3 MeV are selected inside the full neutron trigger gate of 
[16, 1600] us after the muon event. Also 10 C candidates must satisfy a cut in visible energy 
of [2, 4] MeV and occur within [10, 310] s to a preceding ^n-coincidence. A lower energy 
threshold of 2 MeV avoids a contribution of 11 C decays in the data set. The energy dis- 
tribution of the 10 C candidates is constructed from events with At G [10, 50] s. Only 11 Be 
contributes as a cosmogenic contaminant in this parameter selection. Based on the selec- 
tion cuts and the additional requirement of a ^m-coincidence, the contribution of 11 Be is 
estimated to be less than 6 • 10 -3 (d 100t) _1 and this is taken into account as a system- 
atic uncertainty. The spectral shape of uncorrelated background is derived from events at 
At > 310 s. The simultaneous fit is depicted in figure 10. The fraction of 10 C decays 
accompanied by a muon in coincidence with at least one detected neutron capture is esti- 
mated via a test sample of 10 C candidates. 10 C candidates are selected by At £ [10, 310] s 
and a lateral distance of 1 m to a parent muon after removal of the neutron requirement. 
The numbers of 10 C decays in the subset which satisfy the neutron requirements (subset 
A), as well as in the complementary subset (subset B), are derived by time profile fits and 
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Figure 9. Simultaneous fit of the cosmogenic isotope u Be in visible energy deposition (top panel) 
and decay time relative to preceding muons (bottom panel). The fit returns only an upper limit for 
the isotope n Be. The goodness of the simultaneous fit is x 2 / n df = 30/61. We note, that the low 
value of the reduced \ 2 is connected to the low statistics of the data set fit in the unbinned maximum 
likelihood fit. 



then compared. The 10 C tagging efficiency due to the neutron requirement is evaluated 
to e n ( w C) = (92.5]t2Q.o) The broad uncertainty range is associated with the determi- 
nation of 10 C decays in subset B, which is dominated by physically uncorrelated matches. 
Based on the Monte Carlo simulation, we expect e( 10 C) = (79.0 ± 0.5) % of all 10 C decays 
within the selected energy range. By considering these corrections, the 10 C rate and yield 
are determined by the simultaneous fit to i?( 10 C) = (0.52 ± 0.07 sta tta06syst) (dlOOt)^ 1 and 
Y( 10 C) = (1.79 ± 0.25 stat i°;| syst ) • lO" 6 /(p ■ (g/cm 2 )). 

5.7 n C 

Like 10 C, neutron emission is expected in the muon-induced production of cosmogenic n C 
which is a /3 + -emitter, with r = 29.4 min , Q = 1.98 MeV. Therefore, an analogous three-fold 
coincidence between a muon, subsequent neutron capture(s), and 11 C candidates is applied 
in the rate determination. Candidates of n C decays are selected within the energy range of 
[1,2] MeV and the time gate of [0.1, 3.6] h with respect to a preceding muon- neutron coin- 
cidence. The events for the energy distribution are selected in the same time interval. As 
a result of the long lifetime of n C and an average run duration of ~6h in Borexino, effects 
of run boundaries on the time profile are not negligible. To avoid a distortion of the time 
profile, n C decays within the first 3.6 h after run start are not considered in the analysis. 
This restriction reduces the data set to a live time of 188 d. Events within 2h and 4 m with 
respect to any neutron capture vertex of a /m-coincidence are vetoed to obtain the formation 
of the background spectral shape. The remaining events are used to derive a spectrum con- 
taining only ^n-uncorrelated background sources, i.e. non-cosmogenic background, and n C 
production without the detection of a subsequent neutron capture. The latter contribution 
is the result of the limited neutron detection efficiency, in case of saturated detector electron- 
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ics and so-called invisible channels. Invisible channels denote all muon-induced production 
processes, yielding n C with no free neutron emission in the final state [9, 16]. 

The result of a simultaneous fit in energy and time is shown in figure 11. The fraction of 
11 C decays correlated with a muon and at least one detected neutron is estimated in the same 
manner as in the 10 C analysis described in section 5.6. The selection of the subset of 11 C 
candidates is chosen within At S [0.1, 3.6] h and a lateral distance of 1 m to a preceding \in- 
coincidence. The obtained efficiency for the neutron requirement is e ra ( 11 C) = (86.8 ±6.9) %. 
The cosmogenic production rate is found to be R( n C) = (25.8 ± 1.3 stat ± 3.2 syst ) (dlOOt) -1 
with a corresponding yield of Y( n C) = (8.86 ± 0.45 sta t ± l-10 syst ) • 1(T 5 /{jj, ■ (g/cm 2 )) after 
correcting for the fraction of e( 11 C) = (92.2 ± 0.4) % of all decays which deposit a visible 
energy in the selected parameter space. The rate of 11 C decays detected in coincidence 
with cosmic muons and associated neutrons is in good agreement with results obtained from 
spectral fits without the coincidence requirement. Modeling all signal components in the 
energy range of [270, 1600] keV, the n C rate was found to be R( U C) = (28.5 ± 0.2 stat ± 
0.7 sys t) (d 100t) _1 in the precision measurement of the solar 7 Be neutrino interaction rate in 
Borexino [2]. 

6 Monte Carlo simulations of Cosmogenic Neutrons and Radioisotopes 

The production of cosmogenic neutrons and radioisotopes in Borexino was studied by using 
the Geant4 [17, 18] and Fluka [19, 20] simulation packages which are both commonly used 
to simulate deep underground, low background experiments. Their predictions are compared 
with our experimental data. 
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Figure 11. Simultaneous fit of the cosmogenic isotope 11 C in visible energy deposition (top panel) 
and decay time relative to preceding muons (bottom panel). The goodness of the simultaneous fit is 
X 2 /ndf= 189/168. 



6.1 Simulation procedure 

Initially, a careful description of the muon-induced radiation field at LNGS was prepared 
using Fluka. For this simulation a model of Hall C is surrounded by a 700 cm thick shell of 
Gran Sasso rock [21] which is found sufficient to allow a full shower development. The setup 
is subjected to residual cosmogenic muons taking into account the muon angular distribution 
and the muon differential energy spectrum as a function of the slant depth and muon event 
multiplicity as measured by the MACRO experiment [12, 22]. The adopted muon charge 
ratio is = N^+fN^- ~ 1.38 as measured by the OPERA experiment [23]. Cosmogenic 
muons and muon-induced secondaries emerging into Hall C are followed, and all particles 
reaching the Borexino water tank are recorded. Details of this simulation are given in [24]. 

A fraction of 1.5 % of the cosmogenic events with multiple muons crossing the Borexino 
detector simultaneously is found from simulation. Further, 12 % of single muon events are 
actually caused by single muons which belong to muon bundles. 

The cosmogenic events which were recorded at the outside of the water tank are then 
used as a source for both Geant4 and Fluka to simulate the production of cosmogenic neu- 
trons and radioisotopes inside the Borexino detector setup. The yields are extracted directly 
by recording neutron captures and residual isotope production, rather than simulating the 
full detector response. The Monte Carlo simulation predictions are compared to the data 
analysis results of the previous sections that include corrections for detection efficiency. 

For the comparison between predictions and experimental data, the following observ- 
ables are considered: 1) the rates for neutron captures and cosmogenic isotopes production, 
2) the neutron capture time, 3) the neutron capture multiplicity for individual muon events, 
and 4) the lateral distance between the neutron capture and the parent muon track. We also 
compare the rate of muon events for which one or more neutron captures are recorded. 
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6.1.1 Geant4 

The version of Geant4 which is used for this study is Geant4-09-06-patch-01 from February 
2013. Detailed description of the Geant4 toolkit is available in [25]. 

Many hadronic models are available in Geant4 [26] , and can be used as function of the 
particle energy. We cover the entire muon energy spectrum at LNGS through a combination 
of different physics inputs. For our study, we use: a) the Quark-Gluon String (QGS) model 
for proton, neutron, pion and kaon interactions with nuclei at kinetic energies above 12 GeV, 
completed with the Precompound model for the evaporation phase of the interaction ; b) 
the Fritiof model (FTF) for the interaction of highly energetic protons, neutrons, pions and 
kaons starting from 4-5 GeV, also completed with the Precompound model ; c) the Bertini 
cascade (BERT) model, which includes intra-nuclear cascade, followed by precompound and 
evaporation phases of the residual nucleus, for proton, neutron, pion and kaon interactions 
with nuclei at kinetic energies below 9.9 GeV ; d) the Binary cascade (BIC) model, a data 
driven intra-nuclear cascade model intended for energies below 5 GeV ; e) the High Precision 
Neutron (HP) model, describing parameterized capture and fission for low-energy neutrons 
(below 20 MeV) ; f) the Low Energy Parameterized (LEP) model for proton, neutron, pion 
and kaon interactions with nuclei at kinetic energies between 9.5 GeV and 25 GeV. 

Some "ready-made" Physics Lists merging different models are available, and we have 
defined four models: Model I (merge of FTF and BIC with HP manually added), Model II 
(merge of FTF, BERT and HP), Model III (merge of QGS, BIC and HP), Model IV (merge 
of QGS, BERT and HP). A summary of the introduced hadronic models is shown in table 3. 
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HP 


Binary 


Bertini 


FTF 


Protons 






0^5 GeV 






4 GeV -> 100 TeV 


Neutrons 


- 


-* 20 MeV 


19.9 MeV -» 5 GeV 






4 GeV -> 100 TeV 


TV 
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4 GeV -> 100 TeV 


K 











-> 5 GeV 


4 GeV -> 100 TeV 


Model II 




HP 


Bertini 






FTF 


Protons 






0^5 GeV 






4 GeV -> 100 TeV 


Neutrons 


- 


-> 20 MeV 


19.9 MeV -» 5 GeV 






4 GeV -> 100 TeV 


7T, K 






-> 5 GeV 






4 GeV -> 100 TeV 


Model III 




HP 


Binary 




LEP 


QGS 


Protons 






-> 9.9 GeV 


9.5 


-> 25 GeV 


12 GeV -> 100 TeV 


Neutrons 


- 


-> 20 MeV 


19.9 MeV -> 9.9 GeV 


9.5 


-> 25 GeV 


12 GeV -> 100 TeV 


7T, K 






-> 9.9 GeV 


9.5 


->• 25 GeV 


12 GeV -> 100 TeV 


Model IV 




HP 


Bertini 




LEP 


QGS 


Protons 






-> 9.9 GeV 


9.5 


->• 25 GeV 


12 GeV -> 100 TeV 


Neutrons 


- 


-> 20 MeV 


19.9 MeV -> 9.9 GeV 


9.5 


-> 25 GeV 


12 GeV -> 100 TeV 


7T, K 






-> 9.9 GeV 


9.5 


-> 25 GeV 


12 GeV -> 100 TeV 



Table 3. Summary of the Hadronic Models used in Geant4. 

Each model has been tested with and without the activation of the Light Ion (LI) 
Physics List (which defines the light ions likely to be produced by the hadronic interac- 
tions, such as deuterons, tritons, 3 He, a-particles and generic ions). In addition, the fol- 
lowing Physics Lists are included for each model: the Electromagnetic Processes for Lep- 
tons (G4eMultipleScattering, G4eIonisation, G4eBremsstrahlung, ...), the Nuclear Decay 
Processes (G4Decay, G4RadioactiveDecay) and the standard Elastic Scattering for hadrons 
(G4HadronElasticPhysics) . 
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6.1.2 Fluka 



Fluka is a fully integrated particle physics Monte Carlo simulation package based predom- 
inantly on original and well-tested microscopic models. The models are benchmarked and 
optimized by comparing to experimental data at the single interaction level. The physics 
models in Fluka are fully integrated into the code and no modifications or adjustments are 
available at the user level. A list of benchmark results relevant to the simulation of deep 
underground cosmogenic backgrounds are described in [24]. 

Details of the physics models implemented in Fluka with focus on hadronic interactions 
and the Fluka specific nuclear interaction model PEANUT can be found in [27-30]. A 
validation of the Fluka Monte Carlo code for predicting induced radioactivity is given for 
instance in [31]. 

In general, the simulation was performed using the Fluka default setting PRECI- 
SION). In addition, photonuclear interactions were enabled through the Fluka option 
PHOTONUC and a more detailed treatment of nuclear de-excitation was requested with 
the EVAPORAT(ion) and COALESCE(nce) options. These enable the evaporation of heavy 
fragments (A > 1) and the emission of energetic light fragments, respectively. The treatment 
of nucleus-nucleus interaction was turned on for all energies via the option IONTRANS and 
radioactive decays were activated through the option RADDECAY. 

The version of Fluka used for the present study is FLUKA2011.2, released in November 
2011. Further information about the implemented physics models is available through the 
Fluka manual and additional documentation and lecture notes located at the official Fluka 
website [32]. 

6.2 Simulation results 

We present the predictions obtained with Geant4 and Fluka regarding different physics 
observables. For Geant4, we concluded that the best match to data is given by Model III 
or IV, depending on the observable under study. Model I produced 15 % less neutrons than 
the other models and Model II shows no relevant difference with respect to model IV. For 
readability we choose to present here only Model III and IV, however complete results for all 
four models are available as supplementary material of this article. 

6.2.1 Cosmogenic Radioisotopes 

Production yields for cosmogenic isotopes and neutrons are summarized in table 4. The 
yields measured by Borexino are compared to Geant4 and Fluka predictions as well as 
measured yields from the KamLAND experiment [9] given in the rightmost column. We 
note that KamLAND has a different number of carbon nuclei per ton of liquid scintillator 
(4.30- 10 28 for KamLAND as opposed to 4.52- 10 28 in case of Borexino) and that the mean 
residual muon energies differ somewhat between the two sites: (283 ± 19) GeV at LNGS 2 
versus (260 ±8) GeV at the Kamioka mine [9], which should lead to a difference in the 
observed production yields. 

2 The residual mean muon energy is based on the MACRO measurement for single and double muon events 
reported in [33]. Details of the procedure are given in [24]. 
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Geant4 


Geant4 


Fluka 


Borexino 


KamLAND 




Model III 


Model IV 












- (E M ) = 289 ± 19 GeV - 




(E M ) = 260 ± 8 GeV 


Isotopes 




Yield [10" 7 (^! 


i/cm 2 V 1 ] 




12 N 


1.11 ±0.13 


o.U ± U.z 


u.o ± u.z 


< 1.1 


1.8 ±0.4 


12 B 


30.1 ±0.7 


29.7 ±0.7 


28.8 ±1.9 


56 ±3 


42.9 ±3.3 


8 He 


<0.04 


0.18 ±0.05 


0.30 ±0.15 


< 1.5 


0.7±0.4 


9 Li 


0.6 ±0.1 


1.68±0.16 


3.1±0.4 


2.9 ±0.3 


2.2 ±0.2 


8 B 


0.52 ±0.09 


1.44 ±0.15 


6.6 ±0.6 


14 ±6 


8.4 ±2.4 


6 He 


18.5 ±0.5 


8.9 ±0.4 


17.3 ±1.1 


38 ±15 


not reported 


8 Li 


27.7 ±0.7 


7.8 ±0.4 


28.8 ±1.0 


7±7 


12.2 ±2.6 


9 C 


0.16 ±0.05 


0.99 ±0.13 


0.91 ±0.10 


< 16 


3.0±1.2 


n Be 


0.24 ±0.06 


0.45 ±0.09 


0.59 ±0.12 


<7.0 


1.1±0.2 


io c 


15.0 ±0.5 


41.1 ±0.8 


14.1 ±0.7 


18 ±5 


16.5 ±1.9 


n C 


315 ±2 


415 ±3 


467 ± 23 


886 ±115 


866 ±153 


Neutrons 




Yield [10- 4 (^j 


y/cm 2 )^ 1 ] 






3.01 ±0.05 


2.99 ±0.03 


2.46 ±0.12 


3.10±0.11 


2.79 ±0.31 



Table 4. Predicted yields for cosmogenic products obtained from Geant4 (Model III and IV) and 
Fluka are compared to data from Borexino . Also shown are results from the KamLAND experiment 
[9] . Note that the production yields depend on the number of carbon atoms per weight and the muon 
energy spectrum. Thus, a 10-20% difference between KamLAND and Borexino results is expected. 

6.2.2 Cosmogenic neutrons 

Neutron capture time. The simulated neutron capture time of the Borexino scintillator 
from Geant4 and Fluka are (275.8 ± 0.9) us 3 and (253.4 ± 0.6) us, respectively. This is 
to be compared to the measured capture time of (259. 7± 1.3 s t a t ± 2.0 sys t) us. The neutron 
capture time was also measured in Borexino using an Am-Be neutron source [11] which yields 
(254.5 ± 1.8) us. The experimental disagreement with the value measured from cosmogenic 
neutrons could be explained by a fraction of neutrons which are captured on iron in the 
source capsule. This was also observed by KamLAND [9]. 

Neutron production yield. In table 4, the neutron production yield is reported. The 
observed neutron production deficit of the Fluka simulation was studied in [24]. The main 
cause of the deficit was found to be the low cosmogenic production rate predicted for 11 C 
(table 4). At the LNGS depth, the production of 11 C in liquid scintillator is followed by a 
neutron emission in 95% of all cases as was shown by [16]. Since the measured 11 C rate is 
almost 30% of the neutron production rate, and the 11 C rate given by Fluka is roughly 
50 % of the measured value, a reduction of the number of predicted cosmogenic neutrons 
in the order of 15% is expected. The origin of the low n C production rate in Fluka is 
addressed by improvements to the Fermi break-up model [34, 35] which will be available 
with the next Fluka release. The impact of the improved model for the 11 C production 
in liquid scintillator at LNGS energies is currently under investigation. In addition, Fluka 
predicts the production of energetic deuterons (-©kin > 50 MeV) inside the liquid scintillator 

3 The out-dated Geant4 version 4.9.2.p02 returns (254.9 ± 0.6) us and is thus in agreement with the mea- 
sured value. No explanation has been found for the discrepancy between the different Geant4 versions. 
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lateral distance n-p, [m] 

Figure 12. Lateral distance between neutron capture points and the parent muon track: comparison 
of Borexino data to predictions obtained with Geant4 - Model III and IV and Fluka . 

with a rate approximately equal to 4 % of the cosmogenic neutron production rate. Energetic 
deuterons can re-interact with the scintillator leading to a break up (and re-capture) of the 
neutrons. However, these processes are not described by Fluka as no interaction model is yet 
implemented for these deuterons. This further reduces the predicted number of cosmogenic 
neutrons. 

The neutron yield of the Geant4 simulation is in good agreement with the data, but 
the 11 C rate is also ~50% of the measured value. As 11 C is usually produced together with 
a neutron, this indicates that the yields of other neutron production channels are too high 
in Geant4. 

The neutron production yield measured by the KamLAND experiment is (2.79 ± 0.31) • 
10 _4 n//i • (g/cm 2 ) [9]. The latest results from the LVD experiment, which is also located at 
LNGS, indicate a neutron yield of (2.9 ± 0.6) • 10~ 4 n/^ • (g/cm 2 ) [36]. 

Muon-neutron lateral distance distribution. The lateral distance distribution of the 
neutron capture location from the parent muon track is shown in figure 12 and compared 
to predictions from simulation. As described in section 4.4, a radial cut of less than 4 m 
for neutrons and an impact parameter cut less than 4 m for muons were applied in the sim- 
ulations. Moreover, the muon track reconstruction uncertainties were applied a posteriori 
to the simulated distribution. They dominate the shape at small distances. The simulated 
distributions were scaled to match the number of measured neutron captures in order to 
compare the shape of the histograms. The experimental distribution is well reproduced by 
both Geant4 - Model IV and Fluka out to large distances. Geant4 - Model III instead 
reproduces data less accurately. The small differences present at the far range can be at- 
tributed to the additional cuts described in section 4.4, which have not been implemented in 
the simulation. They mostly suppress mis-reconstructed neutrons following showering muons 
which are expected to dominate the distribution at large distances from the track. 
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Figure 13. Neutron capture multiplicity: comparison of Borexino data to predictions obtained with 
Geant4 - Model IV and Fluka . Simulated curves were modified to account for the limited neutron 
detection window [30— 1590] us after the muon. Geant4 - Model III curve does not differ appreciably 
from Model IV curve and is not shown. 

Neutron capture multiplicity. The neutron capture multiplicity distribution is shown 
in figure 13, for the the experimental spectrum as obtained from SYS1 (section 4.2). The ex- 
perimental distribution is compared to predictions by Geant4 - Model IV and Fluka which 
were both scaled to match the live time of the experimental data set. With the exception of 
events with low multiplicities, good overall agreement between data and simulation is found 
for both packages out to multiplicities of hundreds of neutron captures. The shape of the 
multiplicity distribution is somewhat distorted towards small multiplicities for both Monte 
Carlo simulations and data with respect to the true physical distribution, because neutrons 
are only detected if they are captured between 30 and 1590 us after the muon trigger. 

Rate of neutron-producing muons. The rate of cosmic muons crossing Borexino which 
produce at least one neutron is not well reproduced by both simulation packages. The 
measured rate is 67 ± 1 events per day, while Fluka returns 41 ± 3 per day and Geant4 
returns 42.5 ±0.2 (Model III) and 44.6 ±0.2 (Model IV). This discrepancy is also apparent 
from the multiplicity plot in figure 13 and seems to be associated with low multiplicity events. 
No explanation has been found for this discrepancy. 

7 Conclusions 

The Borexino detector offers a unique opportunity to study cosmic backgrounds at a depth of 
3800 m w.e. at the Gran Sasso underground laboratories. The results are not only essential to 
low-energy neutrino analyses, but are also of substantial interest for direct dark matter and 
Ovf3f3 searches at underground facilities. Based on thermal neutron captures in the scintillator 
target of Borexino, a spallation neutron yield of Y n = (3.10 ±0.11) • 10~ 4 n/(^- (g/cm 2 )) was 
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determined. The lateral distance profile was measured based on the reconstructed parent 
muon tracks and neutron capture vertices. An average lateral distance of A = (81.5 ±2.7) cm 
was found. The data results on neutron yield, multiplicity and lateral distributions were 
compared to Monte Carlo simulation predictions by the Fluka and Geant4 framework and 
are largely compatible. The simulated neutron yield of Fluka shows a deficit of ~20 %, while 
the result of the Geant4 simulation is in good agreement with the measured value. However, 
both simulations should be increased as a result on an underprediction of 11 C production. 

The production rates of several cosmogenic radioisotopes in the scintillator were de- 
termined based on a simultaneous fit to energy and decay time distributions. Results of a 
corresponding analysis performed by the KamLAND collaboration for the Kamioka under- 
ground laboratory [9] are similar to our findings. Moreover, Borexino rates were compared to 
predictions by Fluka and Geant4: While there is good agreement within their uncertainties 
for most isotopes, some cases ( 12 B, n C, 8 Li for both codes and 8 B, 9 Li for Geant4 only) 
show a significant deviation between data and Monte Carlo simulation predictions. 
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